More advanced matters
Note. We have a working script prepared for you to put your answers in, called “Working Script Practical 2.Rmd.
This is the second practical, we will dive into more advanced matters: model selection, subject specific parameters, model diagnostics, and model fit. Note that we cover a lot of ground here, most likely too much to finish within this practical. All parts are modular, so you are invited to work on what interests you most. We continue with the same dataset: the open access dataset from Rowland and Wenzel (2020), using the subset of the affective experiences happy, relaxed, anxious, and depressed, measured six times a day for 40 days in total.
The answers to the practicals are available in the
*_answers.html files.
In all practicals in this workshop, we make use of the
mHMMbayes package, version 1.1.0. You can load the package
like
library(mHMMbayes)In this practical, we will also use the following packages:
# install.packages("coda")
library(ggplot2) # for plotting
library(ggmHMM) # plotting mHMM objects
library(coda) # for assessing convergenceWe will also, again, make use of our helper functions
source('./helpers.R')Subject specific parameters
In this section, we obtain and inspect the subject specific parameters.
Emission distribution
We first need to obtain the subject specific emission distribution
parameters, and save the subject specific parameters in an object. When
focusing on the emission distribution parameters, the function
obtain_emiss() can be used to obtain the subject specific
emission distribution parameters, by setting the input argument
level to "subject". Note that the default of
level is "group", so when left unspecified,
the function obtain_emiss() will return the group level
emission distribution parameters.
Exercise 1 Obtain the subject specific emission
distribution parameters using the function obtain_emiss()
and set the input argument level to "subject".
Save the output in the object emiss_subject. Inspect the
output.
# obtaining and saving the subject specific transition probabilities
emiss_subject <- obtain_emiss(out_2st_emotion, level = "subject")Again, to develop an intuition on how much the emission distribution
parameters (e.g., the variable and state specific means) vary over
subjects, it can be helpful to plot the subject specific emission
distribution parameters, for example by adding the subject specific
parameters as dots to a group level barplot of the emission means. This
can be done using the ggmHMM function
plot_emiss().
Exercise 2 Plot (a subset of) the subject specific
emission distribution means, using the function
plot_emiss() from the ggmHMM package and
investigate how much the means differ over subjects. Hint: both
functions have the argument subject, which can be used to
specify a vector of subjects to plot random effects for.
plot_emiss(out_2st_emotion, subject_effects = TRUE, line = FALSE, subject = 1:50) # using ggmHMM function
In addition, it is important to check that the states represent the same construct over subjects. Working with our current dataset, that would mean that for each subject, state 1 represents a ‘good mood’ state, and state 2 represents a ‘bad mood’ state. For example, for the indicator happy, that would translate to the mean of state 1 being higher compared to state 2. Visually, this can be checked by connecting the subject specific means over the states within one indicator by lines.
Exercise 3 Plot (a subset of) the subject specific
emission distribution means connecting the subject specific means, and
visually investigate whether the state 1 is a relatively good state and
state 2 a relatively bad state over subjects. Hint: this can be done by
specifying line = TRUE in plot_emiss().
## plot first 50 subjects to aid interpretability
plot_emiss(out_2st_emotion, subject_effects = TRUE, line = TRUE, subject = 1:50) # using ggmHMM function
Transition probabilities gamma
First, we need to obtain the subject specific parameters for the
emissions and the transition probabilities. When focusing on the
transition probabilities gamma, the function obtain_gamma()
can also be used to obtain the subject specific transition
probabilities, by setting the input argument level to
"subject". Note that the default of level is
"group", so when left unspecified, the function
obtain_gamma() will return the group level transition
probabilities.
Exercise 4 Obtain the subject specific transition
probabilities using the function obtain_gamma() and set the
input argument level to "subject". Save the
output in the object gamma_subject. Inspect the output.
# obtaining and saving the subject specific transition probabilities
gamma_subject <- obtain_gamma(out_2st_emotion, level = "subject")
# inspecting the output, restring the returned output to the first 10 subjects
gamma_subject[1:10]## $`Subject 1`
## To state 1 To state 2
## From state 1 0.598 0.402
## From state 2 0.368 0.632
##
## $`Subject 2`
## To state 1 To state 2
## From state 1 0.737 0.263
## From state 2 0.147 0.853
##
## $`Subject 3`
## To state 1 To state 2
## From state 1 0.96 0.04
## From state 2 0.40 0.60
##
## $`Subject 4`
## To state 1 To state 2
## From state 1 0.716 0.284
## From state 2 0.282 0.718
##
## $`Subject 5`
## To state 1 To state 2
## From state 1 0.922 0.078
## From state 2 0.493 0.507
##
## $`Subject 6`
## To state 1 To state 2
## From state 1 0.767 0.233
## From state 2 0.361 0.639
##
## $`Subject 7`
## To state 1 To state 2
## From state 1 0.438 0.562
## From state 2 0.141 0.859
##
## $`Subject 8`
## To state 1 To state 2
## From state 1 0.948 0.052
## From state 2 0.232 0.768
##
## $`Subject 9`
## To state 1 To state 2
## From state 1 0.615 0.385
## From state 2 0.215 0.785
##
## $`Subject 10`
## To state 1 To state 2
## From state 1 0.839 0.161
## From state 2 0.479 0.521Note that it can be quite insightful to compare the subject specific transition probabilities to the plotted state sequences (see practical 1). Subjects with large self-transitions (i.e., probabilities that denote a transition to the same state as the current state, provided on the diagonal of the transition probability matrix) will probably remain for a long time within a state, while low self-transition probabilities most likely will translate to state sequences where a state is only visited for a short consecutive sequence.
To develop an intuition on how much the transition probabilities vary over subjects, it can be helpful to plot the subject specific state transition probabilities, for example in a heat map with one row per subject, or in stacked dots with a dot for each subject specific transition probability, and columns for each of the possible state transitions.
Exercise 5 Plot (a subset of) the subject specific transition probabilities, and investigate how much the transition probabilities differ over subjects. This can be done in 2 ways:
- With the helper function
PlotTrans(), a barplot may be plotted showing the group level effects, with points for the subject-specific effects (or a subset thereof). It uses the following arguments:- model: a model of class
mHMM - subject: an optional vector of subject to plot subject-specific effects for. If not specified, all subject-specific effects are plotted.
- model: a model of class
- With the function
plot_gamma()from the ggmHMM package, a heatmap may be plotted for each subject (or a subset) separately. This can be done by- setting the input argument
subject_effectstoTRUE, and - setting the input argument
subjectto a vector of subject IDs (e.g.,1:10), or a single number do specify the subject. - specify
ncol_facetto setup the layout
- setting the input argument
# Transition probabilities of first subject using helper function
PlotTrans(out_2st_emotion, subject = 1:50) # plot
# Transition probabilities of first 10 subjects using ggmHMM function
plot_gamma(out_2st_emotion, level = "subject", subject = 1:10, ncol_facet = 5) # plot with facets
Exercise 6 Interpret the transition probabilities. Is there much heterogeneity between persons? Are subjects generally more likely to stay in the same state, or switch states? What does this mean, if we take the substantive interpretation of the states into account?
Note that for a two-state model, the information given by the transition probabilities is limited; The diagonals inform us about the probability of remaining in a state (and therefore the stability/duration of states), but the off-diagonals do not provide us with unique information; if a subject goes out of the state, they only have one left to go to. A three-state model may have more interesting information to provide. It could tell us, for example, whether a subject is more likely to immediately switch between a ‘good’ and a ‘bad’ state, or whether they generally move through a ‘medium’ state first.
Model selection
Fit models with a plausible number of states
One of the challenges when using hidden Markov models is the number
of states to choose. Here, we will limit ourselves to models with 2, 3,
and 4 states. As checking model fit is discussed at a later point in the
practical, we will investigate the different models based on the state
composition and the AIC. However, it should be noted that convergence
and model fit are important characteristics that should be considered
when choosing a model to interpret. The 2-state model was fitted in the
first practical, and saved in the object
out_2st_emotion.
Exercise 7 To save time we have already pre-run a
3-state and 4-state model. Use the function readRDS() to
load the fitted models from the files out_3st_emotion.rds
and out_4st_emotion.rds, and save them in the objects
out_3st_emotion and out_4st_emotion. The
settings we used are shown below.
out_3st_emotion <- readRDS("./data/out_3st_emotion.rds")
out_4st_emotion <- readRDS("./data/out_4st_emotion.rds")####################
# 3 state model
####################
## general model properties
m <- 3
n_dep <- 4
## starting values for gamma
start_gamma_3st <- matrix(c(0.8, 0.1, 0.1,
0.1, 0.8, 0.1,
0.1, 0.1, 0.8), byrow = TRUE, ncol = m, nrow = m)
## starting values for the emission distribution
start_emiss_3st <- list(matrix(c(80, 10, # happy
50, 10,
20, 10), byrow = TRUE, ncol = 2, nrow = m),
matrix(c(80, 10, # relaxed
50, 10,
20, 10), byrow = TRUE, ncol = 2, nrow = m),
matrix(c(10, 5, # anxious
30, 10,
50, 15), byrow = TRUE, ncol = 2, nrow = m),
matrix(c(10, 5, # depressed
30, 10,
50, 15), byrow = TRUE, ncol = 2, nrow = m))
## specifying weakly informative prior for continuous emission distributions
emotion_prior_emiss_3st <- prior_emiss_cont(
gen = list(m = m, n_dep = n_dep),
emiss_mu0 = list(matrix(c(80, 50, 20), nrow = 1), # happy
matrix(c(80, 50, 20), nrow = 1), # relaxed
matrix(c(10, 30, 50), nrow = 1), # anxious
matrix(c(10, 30, 50), nrow = 1)), # depressed
emiss_K0 = rep(list(1), n_dep),
emiss_V = rep(list(rep(5^2, m)), n_dep),
emiss_nu = rep(list(1), n_dep),
emiss_a0 = rep(list(rep(1.5, m)), n_dep),
emiss_b0 = rep(list(rep(20, m)), n_dep),
)####################
# 4 state model
####################
## general model properties
m <- 4
## starting values for gamma
start_gamma_4st <- matrix(c(0.7, 0.1, 0.1, 0.1,
0.1, 0.7, 0.1, 0.1,
0.1, 0.1, 0.7, 0.1,
0.1, 0.1, 0.1, 0.7), byrow = TRUE, ncol = m, nrow = m)
## starting values for the emission distribution
start_emiss_4st <- list(matrix(c(80, 10, # happy
60, 10,
40, 10,
20, 10), byrow = TRUE, ncol = 2, nrow = m),
matrix(c(80, 10, # relaxed
60, 10,
40, 10,
20, 10), byrow = TRUE, ncol = 2, nrow = m),
matrix(c(10, 5, # anxious
20, 10,
40, 10,
60, 10), byrow = TRUE, ncol = 2, nrow = m),
matrix(c(10, 5, # depressed
20, 10,
40, 10,
60, 10), byrow = TRUE, ncol = 2, nrow = m))
## specifying weakly informative prior for continuous emission distributions
emotion_prior_emiss_4st <- prior_emiss_cont(
gen = list(m = m, n_dep = n_dep),
emiss_mu0 = list(matrix(c(80, 60, 40, 20), nrow = 1), # happy
matrix(c(80, 60, 40, 20), nrow = 1), # relaxed
matrix(c(10, 20, 40, 60), nrow = 1), # anxious
matrix(c(10, 20, 40, 60), nrow = 1)), # depressed
emiss_K0 = rep(list(1), n_dep),
emiss_V = rep(list(rep(5^2, m)), n_dep),
emiss_nu = rep(list(1), n_dep),
emiss_a0 = rep(list(rep(1.5, m)), n_dep),
emiss_b0 = rep(list(rep(20, m)), n_dep),
)# Fitting the 3 state model
m <- 3
set.seed(42)
out_3st_emotion <- mHMM(s_data = emotion_mHMM,
data_distr = "continuous",
gen = list(m = m, n_dep = n_dep),
start_val = c(list(start_gamma_3st), start_emiss_3st),
emiss_hyp_prior = emotion_prior_emiss_3st,
mcmc = list(J = 500, burn_in = 200))
# Fitting the 4 state model
m <- 4
set.seed(42)
out_4st_emotion <- mHMM(s_data = emotion_mHMM,
data_distr = "continuous",
gen = list(m = m, n_dep = n_dep),
start_val = c(list(start_gamma_4st), start_emiss_4st),
emiss_hyp_prior = emotion_prior_emiss_4st,
mcmc = list(J = 500, burn_in = 200))Exercise 8 Inspect the general output of the 3- and
4- state models by using the inbuilt print() and
summary() for objects of class mHMM returned
by the function mHMM().
# 3 state model
out_3st_emotion## Number of subjects: 125
##
## 500 iterations used in the MCMC algorithm with a burn in of 200
## Average Log likelihood over all subjects: -2823.466
## Average AIC over all subjects: 5706.932
## Average AICc over all subjects: 5710.502
##
## Number of states used: 3
##
## Number of dependent variables used: 4
##
## Type of dependent variable(s): continuoussummary(out_3st_emotion)## State transition probability matrix
## (at the group level):
##
## To state 1 To state 2 To state 3
## From state 1 0.864 0.087 0.050
## From state 2 0.314 0.575 0.111
## From state 3 0.193 0.110 0.697
##
##
## Emission distribution ( continuous ) for each of the dependent variables
## (at the group level):
##
## $happy
## Mean SD
## State 1 68.106 14.927
## State 2 50.765 22.294
## State 3 40.561 14.148
##
## $relaxed
## Mean SD
## State 1 62.515 20.896
## State 2 45.273 25.916
## State 3 38.207 16.530
##
## $anxious
## Mean SD
## State 1 5.357 3.538
## State 2 27.273 26.663
## State 3 23.466 14.509
##
## $depressed
## Mean SD
## State 1 9.798 7.566
## State 2 36.392 26.328
## State 3 43.209 18.970# 4 state model
out_4st_emotion## Number of subjects: 125
##
## 500 iterations used in the MCMC algorithm with a burn in of 200
## Average Log likelihood over all subjects: -2772.047
## Average AIC over all subjects: 5632.095
## Average AICc over all subjects: 5639.906
##
## Number of states used: 4
##
## Number of dependent variables used: 4
##
## Type of dependent variable(s): continuoussummary(out_4st_emotion)## State transition probability matrix
## (at the group level):
##
## To state 1 To state 2 To state 3 To state 4
## From state 1 0.785 0.145 0.053 0.017
## From state 2 0.234 0.664 0.063 0.038
## From state 3 0.265 0.222 0.449 0.064
## From state 4 0.040 0.179 0.063 0.718
##
##
## Emission distribution ( continuous ) for each of the dependent variables
## (at the group level):
##
## $happy
## Mean SD
## State 1 69.306 14.580
## State 2 55.936 14.678
## State 3 52.343 21.409
## State 4 33.437 13.849
##
## $relaxed
## Mean SD
## State 1 63.676 20.775
## State 2 51.315 19.587
## State 3 45.295 25.826
## State 4 32.520 15.771
##
## $anxious
## Mean SD
## State 1 2.803 2.512
## State 2 10.563 8.595
## State 3 36.242 28.347
## State 4 27.706 20.357
##
## $depressed
## Mean SD
## State 1 6.349 5.129
## State 2 27.410 17.391
## State 3 35.180 26.757
## State 4 52.370 18.392Inspect the state composition of the fitted models
First, we will investigate the composition of the states, and use this to evaluate whether it seems sensible to add extra states. Especially with multivariate data, inspecting the composition of the states is done most easily visually.
Exercise 9 Use the function
plot_emiss() from the ggmHMM package to visualize
the emission distributions of the three-state and four-state model. Do
not plot subject-specific effects. Does it seem to make sense to add
extra state(s)?
plot_emiss(out_3st_emotion, subject_effects = FALSE) + ylim(c(0,100)) # three-state model with helper function
plot_emiss(out_4st_emotion, subject_effects = FALSE) + ylim(c(0,100)) # four-state model with ggmHMM function
# Note: the 4th state differentiates the emotions even more. Next, let’s compare the 2-, 3-, and 4-state model on the model selection criterion AIC.
Exercise 10 Obtain the AIC values of the 2-, 3-, and
4-state model using the print() function on the mHMM
objects out_2st_emotion, out_3st_emotion,
out_4st_emotion, and compare. It can be helpful to
visualize the AIC values over the states to get a feeling for how much
the AIC decreases over the models.
out_2st_emotion## Number of subjects: 125
##
## 500 iterations used in the MCMC algorithm with a burn in of 200
## Average Log likelihood over all subjects: -2870.816
## Average AIC over all subjects: 5777.631
## Average AICc over all subjects: 5778.915
##
## Number of states used: 2
##
## Number of dependent variables used: 4
##
## Type of dependent variable(s): continuousout_3st_emotion## Number of subjects: 125
##
## 500 iterations used in the MCMC algorithm with a burn in of 200
## Average Log likelihood over all subjects: -2823.466
## Average AIC over all subjects: 5706.932
## Average AICc over all subjects: 5710.502
##
## Number of states used: 3
##
## Number of dependent variables used: 4
##
## Type of dependent variable(s): continuousout_4st_emotion## Number of subjects: 125
##
## 500 iterations used in the MCMC algorithm with a burn in of 200
## Average Log likelihood over all subjects: -2772.047
## Average AIC over all subjects: 5632.095
## Average AICc over all subjects: 5639.906
##
## Number of states used: 4
##
## Number of dependent variables used: 4
##
## Type of dependent variable(s): continuousAICs <- data.frame(AIC = c(5777.631, 5706.932, 5632.095),
states = c(2:4))
ggplot(AICs, mapping = aes(x = states, y = AIC)) +
geom_point() +
geom_line() +
ggtitle("AIC over the 2 - 4 state models") +
xlab("Number of states") +
theme_minimal()
As you can see, determining the number of states can be a subjective and tricky endeavor, and depends very much on the data and research question at hand.
Model diagnostics and fit
In this section we will inspect model diagnostics and fit, checking convergence, label switching, using posterior predictive checks (PPCs), and evaluating Pseudo residuals.
Convergence
Within any Bayesian analysis, convergence needs to be assessed in order to rule out output resulting from a local (instead of global) maximum. To check this, multiple chains (now we will use 2, but we recommend using at least 4) need to be run, performing the same analysis but using (slightly) different starting values in each chain. Note that in multilevel hidden Markov models, it is important to still keep sensible starting values when varying these. For example, when modelling continuous or count data, we need to keep the same sensible ordering of \(\bar{\mu}_i\) or \(\bar{\lambda}_i\) given data and process. Values of the (weakly informative) prior distributions should be fixed over the chains.
Exercise 11 Specify new starting values to fit at least one additional chain for the 2 state multilevel HMM model. Do not yet run the model.
Exercise 12 To save time, we will provide you with
an additional chain, run with the starting values below. Use the
function readRDS() to load the fitted model from the file
out_2st_emotion_b.rds, and save it in the object
out_2st_emotion_b.
Note. The priors that we used may be slightly different from the ones you specified in your model in the first practical. For now, we will ignore this for simplicity, but it is important to use the same settings for the prior when applying the model yourself.
out_2st_emotion_b <- readRDS("./data/out_2st_emotion_b.rds")m <- 2
n_dep <- 4
start_gamma_b <- matrix(c(0.9, 0.1,
0.1, 0.9), byrow = TRUE, ncol = m, nrow = m)
start_emiss_b <- list(matrix(c(80, 10, # happy
40, 15), byrow = TRUE, ncol = 2, nrow = m),
matrix(c(80, 10, # relaxed
40, 15), byrow = TRUE, ncol = 2, nrow = m),
matrix(c(5, 5, # anxious
30, 15), byrow = TRUE, ncol = 2, nrow = m),
matrix(c(5, 5, # depressed
30, 15), byrow = TRUE, ncol = 2, nrow = m))
out_2st_emotion_b <- mHMM(s_data = emotion_mHMM,
data_distr = "continuous",
gen = list(m = m, n_dep = n_dep),
start_val = c(list(start_gamma_b), start_emiss_b),
emiss_hyp_prior = emotion_prior_emiss,
mcmc = list(J = 500, burn_in = 200))Evidence for non convergence can be checked in multiple ways. A first way is, using trace plots which visualize the parameter estimates over the iterations of the MCMC sampler. The chains should not show any trend and good mixing. Furthermore, we should have no label switching (or at least very little) in the subject-specific parameters (see the section ‘label switching’ later in the practical)..
The trace plots can be obtained using the ggmHMM function
plot_trace(). For now, the following arguments are
relevant:
- model: a model or list of models of class
mHMM - component: “gamma” for transition probabilities, or “emiss” for emission distributions
- level: “group” for group-level parameters. “subject” for subject-specific parameters.
- vrb: the variable to plot the trace for. For example, “happy” for the happy variable.
- subject: Integer specifying the subject to create trace plots for.
Exercise 13
A. Use plot_trace() to create a trace
plot for the group-level transition probabilities of the 2-state model.
Use both out_2st_emotion and
out_2st_emotion_b. Hint: you leed to put the models in a
list using list() to use them in
plot_trace().
plot_trace(model = list(out_2st_emotion, out_2st_emotion_b),
component = "gamma",
level = "group")
B. In addition to the group level parameters, it is
also important to check convergence at the subject level. Use the
function plot_trace() to create a trace plot for the
subject specific emission means of the first subject, for the variable
“happy”. Use both out_2st_emotion and
out_2st_emotion_b.
plot_trace(model = list(out_2st_emotion, out_2st_emotion_b),
component = "emiss",
level = "subject",
vrb = "happy",
subject = 1) + ylim(c(0,100))
We should also inspect convergence numerically, with the potential scale reduction factor \(\hat{R}\). \(\hat{R}\) evaluates equality of means of the different chains - the degree to which variance (of the means) between chains exceeds what one would expect if the chains were identically distributed. A value of R-hat below 1.2 is used to indicate convergence.
We have provided you with a helper function, f_GR(), to
calculate the group level \(\hat{R}\)
values. The function takes the following arguments:
- model_list: a list of models (i.e. multiple chains)
- digits: number of digits to round to. Default is 2.
- burnin: number of burnin iterations to use. Defaults to the number specified when building the first model.
It returns a list with 2 matrices:
- m_RG_gamma is a matrix with m rows and m-1 columns. This is due to the fact that the transition probabilities are estimated using a multinomial logit model, where only m-1 values are estimated (and the last one is derived from the others).
- RG_emiss is a matrix with n_dep rows and m columns. Therefore, the value on the first row and the second column refers to \(\hat{R}\) for the first dependent variable (e.g. “happy”) and the second state.
Exercise 14 Use the function f_GR() to
calculate the \(\hat{R}\) values for
the group level transition probabilities and emission distribution
means. Did the model reach convergence?
f_GR(model_list = list(out_2st_emotion, out_2st_emotion_b),
digits = 2,
burnin = 200)## $m_RG_gamma
## [,1]
## [1,] 1.01
## [2,] 1.06
##
## $RG_emiss
## [,1] [,2]
## [1,] 1.00 1.01
## [2,] 1.07 1.00
## [3,] 1.01 1.06
## [4,] 1.19 1.14Label Switching
A concept related to convergence is that of label switching. In a multilevel HMM, it does not matter for the model what the ordering of the states is. For example, a model where State 1 is the ‘good’ state (e.g. with high “happy” and low “depressed”), and state 2 the ‘bad’ state, or whether it is the other way around. In the iterative MCMC process, the model may switch the labels of the states over iterations, and therefore the interpretation of the states. This is called label switching. In the mHMM, label switching happens at the subject level, but consequently also affects estimates of the group-level effects. It is important to check to what extent this happens.
To check for label switching, it is helpful to use trace plots for
the subject-level parameters, similar to those from last exercise. It
may, however, be useful here to show the trace plots from multiple
states in a single plot, to see whether the states switch labels. This
can be done with the helper function
plot_label_switching().
plot_label_switching() takes as arguments: * model:
output object from mHMM() * subject: a vector of integers
to select the subjects to create plots for. * vrb: a character vector to
specify the variable to plot.
Exercise 15 Use the function
plot_label_switching() to evaluate label switching for
subject 1 to 5, for the 3-state model. Is there a label switching issue
for these persons? And for subjects 6 to 10?
plot_label_switching(out_3st_emotion, subject = 1:5)
plot_label_switching(out_3st_emotion, subject = 6:10)
Exercise 16 The estimate of subject 9’s state specific means for each variable is calculated by taking the mean (or median) of all sampled values after the specified burnin period. How does the label switching issue affect the estimates of this subject’s state-specific means?
PPCs
Finally, we will assess whether the model recovers the data correctly
with respect to the mean and standard deviation, to reveal potential
model misspecification. It is particularly helpful for group-level model
evaluation. First, we need to generate simulated data sets based on
obtained parameter estimates. The function sim_mHMM() can
be used to simulate data using a multilevel hidden Markov model. Note
that for each simulated dataset, new subject-specific parameters are
simulated, but the group-level parameters remain fixed.
The function takes the following arguments (for more detailed
information, use help(sim_mHMM()):
- n_t: length of observed sequence for each subject
- n: sumber of subjects
- data_distr: string indicating the type of data to simulate. “continuous” in our case.
- gen: list containing the elemens: m (number of hidden states) and n_dep (number of dependent variables)
- gamma: group-level transition probability matrix. It should be an unnamed matrix with m rows and m columns.
- emiss_distr: group-level emission distributions. Specified as a list of matrices, one for each dependent variable. Each matrix should have 2 columns (mean and sd) and m rows (number of states).
- var_gamma: variance of transition probabilities.
- var_emiss: variance of emission probabilities.
Note. In the data, there is some observed missingness, including the night gap. For now, we simulate data without the night gap and ignore the missingness. Although it will likely make little difference, it may be sensible to insert missing values in places where the data are missing when performing PPCs on your own data.
Exercise 17 Use the starter script below to simulate 100 datasets based on the model parameters of the 2-state multilevel HMM.
n_sim <- 100 # for now, we will use 100 simulated datasets, but feel free to change.
sim_datasets <- vector("list", n_sim)
m <- ... # specify the number of states
n_dep <- ... # specify the number of dependent variables used in mHMM()
# calculate variance in emission distribution
sapply(out_2st_emotion$emiss_varmu_bar, apply, 2, mean, na.rm = TRUE)
for(i in 1:n_sim){
sim_datasets[[i]] <- sim_mHMM(n_t = 240 , # specify the number of beeps per subject
n = 125, # specify the number of subjects
data_distr = "continuous",
gen = list(m = m, n_dep = n_dep),
gamma = ... , # specify the obtained transition prob matrix
emiss_distr = ... , # specify the emission distribution, see help file
var_gamma = 0.01,
var_emiss = c(300, 300, 200, 200)
)
}n_sim <- 100
sim_datasets <- vector("list", n_sim)
m <- 2
n_dep <- 4
# calculate variance in emission distribution
sapply((out_2st_emotion$emiss_varmu_bar), apply, 2, mean, na.rm = TRUE)## happy relaxed anxious depressed
## varmu_1 323.8858 349.0863 77.74446 174.3764
## varmu_2 324.4923 271.7051 441.86128 314.2728set.seed(42)
for(i in 1:n_sim){
sim_datasets[[i]] <- sim_mHMM(n_t = 240,
n = 125,
data_distr = "continuous",
gen = list(m = m, n_dep = n_dep),
gamma = unname(obtain_gamma(out_2st_emotion)),
emiss_distr = obtain_emiss(out_2st_emotion),
var_gamma = 0.01,
var_emiss = c(300, 300, 200, 200)
)
}Exercise 18 For each simulated data set, calculate the summary statistics of interest for all four variables (i.e., mean and standard deviation).
means_simdata <- matrix(, nrow = n_sim, ncol = n_dep)
for(i in 1:n_sim){
means_simdata[i,] <- apply(sim_datasets[[i]]$obs[,2:(n_dep + 1)], 2, mean)
}
sd_simdata <- matrix(, nrow = n_sim, ncol = n_dep)
for(i in 1:n_sim){
sd_simdata[i,] <- sqrt(apply(sim_datasets[[i]]$obs[,2:(n_dep + 1)], 2, var))
}Now that we have simulated datasets and obtained the summary characteristics, we will visually compare the characteristics of the simulated data to that of the observed data.
Exercise 19 Plot the summary characteristics (mean and sd) of the simulated datasets in histograms, and add a vertical line for the mean and sd in the observed data to aid the comparison.
vars <- names(obtain_emiss(out_2st_emotion))
par(mfrow = c(2,4))
for(i in 1:n_dep){
hist(means_simdata[,i],
main = vars[i],
xlab = paste("Mean", vars[i])
)
abline(col = "red", lwd = 2, v = mean(emotion_mHMM[,1+i], na.rm = TRUE))
}
# plot histogram sd simulated data and observed in real data
for(i in 1:n_dep){
hist(sd_simdata[,i],
main = vars[i],
xlab = paste("SD", vars[i]),
xlim = c(20,30)
)
abline(col = "red", lwd = 2, v = sqrt(var(emotion_mHMM[,1+i], na.rm = TRUE)))
}
# The means show no evidence of bad model fit, however the standard deviations of the simulated datasets are higher compared to the observed data, especially for happy and relaxed. It would be interesting to see whether a 3-state model would partly solve this.Pseudoresiduals
Where PPCs are most useful for group-level model evaluation,
pseudoresiduals are most useful for model evaluation at the subject
level. These are residuals that are obtained based on the most likely
state sequence for each person (hence the suffix pseudo),
using the Viterbi algorithm (see previous practical).
We can compute the Root-Mean-Squared-Error (RMSE) of the pseudoresiduals, which is a measure of how well the model fits the data, where a lower RMSE indicates a better fit. This is especially useful to compare models.
We have provided you with two helper functions for this.
GetResid() can be used to obtain residuals for a specific
person for a specific variable. The function takes the following
arguments:
- data: the dataset used to build the model
- model: the fitted model
- vrb: string specifying the variable obtain pseudoresiduals for
- subject: integer specifying the subject to obtain pseudoresiduals for
It outputs a list with 4 elements:
- emp: the data of the subject for the variable ‘vrb’
- model: the predicted values of the subject for the variable ‘vrb’
- resid: the pseudoresiduals of the subject for the variable ‘vrb’
- RMSE: the RMSE for the subject for the variable ‘vrb’
Exercise 20
A. Use GetResid() to obtain
pseudoresiduals for subject 1 for all variables. Do this for both the
2-state and 3-state model. Hint: use the function lapply()
to loop over the variables.
pseudoresid_2st <- lapply(c("happy", "relaxed", "anxious", "depressed"), # for variable 1 to 4
function(j) GetResid(data = emotion_mHMM, model = out_2st_emotion, vrb = j, subject = 1))
pseudoresid_3st <- lapply(c("happy", "relaxed", "anxious", "depressed"), # for variable 1 to 4
function(j) GetResid(data = emotion_mHMM, model = out_3st_emotion, vrb = j, subject = 1))B. Obtain the RMSE for the pseudoresiduals of subject 1 for all variables, and compare the RMSE of the 2-state and 3-state model. Which model fits better for this person?
RMSE_2st <- sapply(pseudoresid_2st, function(x) x$RMSE)
RMSE_3st <- sapply(pseudoresid_3st, function(x) x$RMSE)
print(RMSE_2st)## [1] 20.04446 22.32644 24.64004 22.43806print(RMSE_3st)## [1] 14.58216 21.99121 24.38652 20.12457## The 3-state model has slightly better fit for subject 1 as it has the lowest RMSE for all variables.After obtaining the pseudoresiduals, we can also plot them separately for each person and each variable to assess: 1. Whether the pseudoresiduals are (roughly) normally. 2. Whether there is a trend over time. 2. Whether there is any remaining dependency between residuals over time (can also be assessed using the autocorrelation). 3. Whether there are any outliers. 4. Whether variance is (roughly) stable over time.
We have provided you with a helper function,
PlotResid(), to plot the pseudoresiduals. It takes the
following arguments:
- model: output model from
mHMM() - data: a data frame used to build the model.
- vrb: character string specifying the variable to plot pseudoresiduals for
- subject: vector of subjects to plot pseudoresiduals for.
Exercise 21 Use PlotRes() to plot the
pseudoresiduals for subject 1 to 3 to for the variable
happy, for the 3-state model. How do you evaluate the model
fit for these subjects?
PlotRes(out_3st_emotion,
data = emotion_mHMM,
vrb = "happy",
subject = 1:3)